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1. Introduction 

The computation of surfaces with prescribed mean curvature is an important problem 
in numerical relativity with an abundance of applications. We will only give a few 
examples here, which we had in mind while developping this algorithm. 

The first example are marginally trapped surfaces, that is closed spherical surfaces 
E in a Riemannian 3-manifold (M, g), such that iJ ± P = on S. Here H denotes the 
mean curvature of E and P = tr^K is the trace of K along S. The extra tensor field K 
on M represents the second fundamental form of M in spacetime. Apparent horizons 
are outermost marginally trapped surfaces and are used in numerical simulations 
for various purposes like black hole location or excision of black holes in numerical 
simulations inside their apparent horizon. For an introduction and furter references 
see for example PE]- 

In the special, time symmetric case K = 0, i.e. P = 0, these surfaces are minimal 
surfaces, that is critical points of the area functional. When surfaces of minimal area 
are considered that satisfy the additional constraint, that they include a given volume 
with the minimal surface, the minimal surface can be "blown up" to render an evenly 
spaced, geometrically defined foliation of the exterior domain, as it is explained in 
Huisken and Yau |H] and Huisken The surfaces of this foliation have constant 
mean curvature 

H = const . 

This the Euler-Lagrange equation to the constrained problem. It is a quasilinear, 
second order, degenerate elliptic equation for the position of the surface. The size of 
the enclosed volume, as well as the area of the constant mean curvature surfaces is 
then a generic candidate for a geometric radial coordinate. 

Besides this, there are more applications for constant mean curvature foliations. 
One of them is to define a concept for the center of mass of an isolated gravitating 
system Another application is that a geometrically defined foliation also can 
be used to construct geometric gauge conditions for time evolution in the Einstein 
equations as described in [S]. In addition the proof of the Riemannian Penrose 
inequality by Bray |H] uses isoperimetric surfaces, which are special cases of constant 
mean curvature surfaces. This proof establishes the monotonicity of the Hawking 
mass on these surfaces with increasing enclosed volume, given a positive energy 
condition. This energy condition translates into the geometric condition for M to have 
nonncgative scalar curvature. The numerically computed constant mean curvature 
surfaces may therfore be used to detect and measure gravitational fields. 

A considerable number of methods have been exploited to find apparent horizons, 
a current overview of which is given in Thornburg 7^ . Baumgarte and Shapiro [2] also 
review some techniques for locating apparent horizons. Schnetter [H| computed the 
more general "constant expansion surfaces" with H zL P = const. These methods 
can be classified into three different approaches, namely flow like methods using a 
curvature flow to locate the apparent horizon, direct methods using a Newton method 
to solve the elliptic apparent horizon equation, and indirect minimization methods 
that try to minimize an L^-error integral. 

The flow like methods are very robust as they converge for a large set of initial 
data and are able to model the change of topology by using level set methods 
PlEi- However, these methods are rather slow. The more and more popular direct 
elliptic methods are very fast but have a small domain of convergence as is noted by 
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Thornburg [7\. The minimization methods are problematic, since minimizing the Li- 
enor functional corresponds to solving a fourth order PDE, while the actual problem 
is only a second order PDE. This does not only spoil the condition of the problem, 
but may also introduce invalid solutions. All these methods, however, share the fact 
that they use finite differencing. 

In contrary, this article describes a finite element based direct minimization 
method to compute constant mean curvature surfaces 

H = const , 

that inherits it the big domain of convergence from the flow like method, while rivaling 
the speed of the direct elliptic method. The method as it is presented here can not 
be applied without modification to horizon finding in the general non-time symmetric 
case with K ^ 0, but an appropriate modification is outlined in section 1X^1 

We consider R'^ to be equipped with a general metric gij. For the purposes 
of general relativity this metric will be asymptotically fiat, that is in rectangular 
coordinates the metric components satisfy 

g.,,=S,,+0{r-^). 

We also will consider two-dimensional spherical surfaces E C R'^. The induced metric 
on E will be denoted by 7, the outer normal by v and the second fundamental form 
A = Vv. The mean curvature is labeled H = tr-^A. In the following we use Einstein's 
summation convention such that Latin indices range from 1 to 3 whereas Greek indices 
range from 1 to 2. We will frequently use the abbreviation CMC for "constant mean 
curvature" . 

Section El will explain the relationship of constant mean curvature surfaces to 
the isoperimetric problem and give some theoretical results concerning existence and 
properties of constant mean curvature surfaces. The algorithm is explained in section|31 
and some numerical examples are given in section 0J 

2. Theoretical Background 

This paper uses the fact that constant mean curvature surfaces are critical points of 
the isoperimetric problem. The isoperimetric problem is to find a surface enclosing a 
given volume that has minimal area. Formulated precisely, this is 

Definition 2.1 Let M be a Riemannian manifold, then a compact subset D, C M is 
called a solution to the isoperimetric problem if for all fl' C M with Vol{fl') = Vol{n) 
the inequality 

\dn\ < \dn'\ 

holds. 

To treat this problem with tools from the calculus of variations one denotes by a 
variation in AI a smooth map F : M x (—£,£) — > M such that F{-,0) : M ^ M is 
the identity and for all t S (— e, e) the map F{-, t) : M ^ M is a diffeomorphism. We 
collect the following facts from the literature, see eg. fo^' ^ background metric 
or ^21 in the general case. 
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Lemma 2.2 Civen any set d M with smooth boundary S dQ and a variation F 
with normal velocity f = g [ ^|t_Q i ^) on E, the following variation formulas for the 
volume Vol{-) of fl and area | • | o/ E hold. 

Voi{F{n,t))= [ fdfis 

t=0 "'S 



dt 
d_ 

dt 



|F(E,i)| = / Hfd^l^ 



t=o Jt. 

Therefore H can be interpreted as L^-gradient of the area functional, and the constant 
function 1 as the i^-gradient of the volume functional. We have the following 
characterization of volume preserving variations. 

Lemma 2.3 If F is a variation that preserves the volume offld M, then fd^j: ~ 
for E = Oil and f the normal velocity of F on Y,. Conversely, if f is a function 
with fdfiY, = then there exists a volume preserving variation with f as normal 
velocity. 

Introducing the usual Lagrangian, the Euler-Lagrange equation of the isopcrimetric 
problem can be computed. 

Proposition 2.4 If is a smooth solution to the isopcrimetric problem, then is 
bounded by a constant mean curvature surface. 

This only characterizes critical points of the isopcrimetric problem. To give a better 
description the usual concept of stability has to be introduced. 

Definition 2.5 A constant mean curvature surface is called stable if the second 
variation of area in volume preserving directions is nonnegative, and strictly stable if 
it is positive. 

Due to the characterization of volume preserving variations in lemma a sufficient 
condition for stability is the nonnegativity of the Jacobi operator 

J/ = -A/-/(|A|2+Ric(i., v)) 

that is the inequality 



/2(|^|2+Ric(j.,i.)) dfi<Jjdf\' 



for all / e C°°(E) with J / = 0. A constant mean curvature surface S is strictly 
stable, if there exist a > such that 



d/z 



a rdfi< / fjfdfi 

for all / with J f — 0. A strictly stable constant mean curvature surface is an isolated 
local minimum of the isoperimetric problem, in the sense that there is no volume 
preserving variation that does not increase area. 

The algorithm presented was constructed having in mind spatial slices of isolated 
gravitating systems in general relativity. These slices have (in absence of linear 
momentum) the following asymptotic behavior of the metric. 
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Definition 2.6 A strongly asymptotically flat manifold is a Riemannian manifold M 
together with a metric g, such that there is a compact set C C M and a diffeomorphism 
X : M\C ^ \ i?fl'(0) for some R and, such that, in the coordinates given by x, the 
metric g has the form 

/ ui \ ^ 

9i3 = + 2^ j ^'^^ 
with the following decay conditions 

\Q^,\<Cr-^ |9'Qy| < Z = 1,2,3,4 

In this setting Huiskcn and Yau 3 have proved the following 

Theorem 2.7 Let [M, g) he a strongly asymptotically flat manifold with to > 0. Then 
there exists a compact set C <Z M such that on M \ C exists a unique foliation by 
spherical constant mean curvature surfaces, such that 

(i) for growing radius these surfaces approximate Euclidean spheres, 

(ii) the centers of these spheres converge to a point in T{? , 

(Hi) and the respective surfaces are strictly stable with respect to the isoperimetric 
problem. 

For the purposes of this article, the way Huisken and Yau establish the existence of 
constant mean curvature surfaces is very interesting. They use the volume preserving 
mean curvature flow. Solving this flow means finding a map F : "Eq x (0, T) — > for 
an initial surface Eq C M with 
dF 

— = (/i - H):^ for t>0 
dt ' 

F{0) = So 

where h — \Et\^^ J^^ Hdjit with = F[YiQ,t). Huisken and Yau show that in the 
above setting a solution exists for all times, in case the flow is started from a Euclidean 
sphere of radius bigger than some critical radius. For t ^ oo the surfaces St then 
converge to a surface with constant mean curvature. 

In view of the variational perspective, this is the flow to the gradient of the area 
functional projected onto the volume preserving variations, which is a technique for 
solving constrained minimization problems. Huisken and Yau show that this method 
converges. The key issue in this proof is the stability of the CMC-surfaces. That is, 
these surfaces are local minimizers of the isoperimetric problem. 

These facts enable us to use a constrained minimization algorithm to approach the 
numerical computation of constant mean curvature surfaces, since the most efficient 
of these methods rely on the positivity of the second variation. 

3. Description of the Algorithm 

As explained before, Huisken and Yau W have shown that a projected gradient 
flow method for constrained minimization converges in the analytic case. From the 
numerical viewpoint there are much better methods for solving such problems, since 
projected gradient methods tend to converge slowly and do not preserve the constraint 
during iteration. 

The approach taken here is to convert the constrained problem into an 
unconstrained minimization problem, which is then solved by a preconditioned 
conjugate gradient method. 
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3.1. The surface model 

Discrete surfaces are modeled as a triangulated meshes, that is sets of vertices, edges 
and triangles linked together according to the topology of the triangulation. 

To each vertex p of a given triangulation T a variation vector r(p) is attached. 
A vertex can only move into this direction. Associated to T is a discrete set of linear 
finite elements 

iS = < ^ ap0p : ttp e R 
[per 

where (f)p (q) — 5pq for all vertices q of the triangulation and 4>p restricted to a triangle 
is linear. As indicated in the definition, the functions 0p form a basis of S, called the 
nodal basis. A function u G 5 is therefore characterized by its values in the vertices 
of the triangulation. 

Given any u € S the graph of u over T is the triangulation T with the vertices 
p — p + u{p)r{p). Fix a reference triangulation T and model the unknown surface as 
graph of a finite element function over T. 

For the purposes of applying hierarchical basis preconditioning, the implementa- 
tion uses a hierarchical data structure described by Leinen This object oriented 
approach stores the triangles in a tree, the roots being coarse triangles with their 
subtriangles as child nodes whenever the coarse triangle is divided. The edges form 
a similar structure. To link edges and triangles, these two trees have cross references 
according to their topology, such as triangles consisting of edges and edges being part 
of triangles. 

An excellent outline of the method of finite elements and issues of numerical 
integration as needed in the next section is given in jl4| . 

3.2. Discrete Area and Volume Functionals 

To compute the area of a triangulation it is clearly possible to compute the area of 
each triangle and sum up. The computation of the area of a single triangle takes place 
in a standard situation. Define the standard triangle Ts C as the set 

Ts — {{xi,X2) e R^ : Xi,X2 > 0,Xi + X2 < 1 } 

Every triangle T — A(pojPi,P2) C R'^ of the triangulation is diffeomorphic to Ts via 
the linear map 

Ft : Ts ~* T : {xi,X2) ^ {I - xi - X2)po + xipi + X2P2 ■ 
Then the area of T is given by 

AiT)^ j y^det(7,^)dx (1) 

JTs 

where 7q,_3 are the components of the induced metric on T in the coordinates induced 
by Ft- If R^ was equipped with the Euclidean metric, ^Jdetiji^ would be a constant 
on T . Then a quadrature formula to integrate exactly is given by 



A{T) = iy^det (70/3(1/3, 1/3)) 



2 

and this is the integration rule used to define the discrete area 

1 



A(T) = -^/det(7a/3(l/3,l/3)). (2) 
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For this integration rule we obtain that 

\A{T)-A{T)\<coh^A{T), 

where h is the (Euchdean) diameter of T and cq = coig, dg, d^g) depends on the metric 
g, in particular, we have the bounds cq < Cr^ in view of definition 12.61 The above 
formula gives that for the whole triangulation we have 

\A{T)-A{T)\<c,hl,^^A{T) 

where /imax is the maximal diameter of the triangles of T. 

The discrete area A is as differentiable with respect to the values of Up as the 
background metric g. Explicit formulas for the derivatives of A can be computed. 
These formulas contain first derivatives of g in the r(p) directions. 

Consider for example the triangle Tq with vertices (po,Pi,P2) and move pi into 
the r{pi) = ri-direction, which gives a family of triangles = A(po,Pi + £fi,P2)- 
Then look at the following map 

Fe -Ts : {xi,X2) ^ (1 - - X2)po + xi{pi + er{pi)) + 0:2^2 

and compute the e-dcrivative of the discrete area expression in using these 
coordinates. This gives 

d_ 

d^ 



^1^^^) = 717— -(722 .9(^1,^1) -712 3(^2, ri)) 

(711^22,1 + 722Aii,i + 27i2Ai2,i) . 



24A(To 

Here Xa=Pa- Po, Jaf3 = 7q/3(1/3, 1/3) and 



p ■ 

Fo(l/3,l/3) 



Similar formulas can be derived for all vertices. 

To define the discrete volume functional, a fixed reference triangulation T and 
a function u is required. The discrete volume functional is defined as the oriented 
discrete volume enclosed by the shell between T and the graph of u over T . Orientation 
is chosen such that positive u gives positive volume and negative u negative volume, 
that is the r[p) arc interpreted to point outward. 

For the computation of the volume corresponding to a single triangle, define the 
standard prism Ps = Ts y- [0,1]. The prism Ps can be mapped to the volume P 
between a triangle T of T with vertices po,pi,P2 and the corresponding triangle T of 
graphii with vertices Pi = Pi + u{pi)r{pi), i = 0, 1, 2 via 

Gp:Ps -^P 

{xi,X2,X3) 1-^ (1 - 2:3) ((1 - Xi - X2)po + Xipi + X2p2) 

+ ((1 - xi- X2) Po + xipi + X2P2) 
as illustrated in figure 13.21 . The oriented volume of P is then given by 



j 1 dvol = j ^det{gij) dx 

= / \ del{gij o Gp) det{dxGp) dx . 



IPs 

We choose the sign of this expression to match the orientation condition given above. 
The integral in this formula is again replaced by a quadrature formula which is exact 
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X, 



Figure 1. Standard prism mapped to volume corresponding to one triangle. 

in the case of a Euclidean background metric. Such a formula can be constructed 
as product of a quadrature formula for Ts and one for [0, 1]. For T5 we can use the 
center of gravity rule previously used for the area. For [0, 1] an integration rule which 
is exact for polynomials of degree two is sufficient. Such a rule is for example given by 
the two-point Gauss rule with weights 1/2 and evaluation points 1/2 ± l/vT2 which 
is exact for polynomials of degree three. The formula for discrete volume computed 
with a general integration rule for [0, 1] with the N weights ak and evaluation points 
Zk then reads 

N 

V{P) - J2 «fc\/det(5,, (Gp(l/3, 1/3, Zk)) det(d,Gp)(l/3, 1/3, Zk) (3) 

k=l 

There is no problem to evaluate these terms, but more elegantly det(d2:G'p)(l/3, 1/3, xs) 
can be written as a polynomial in Up of degree up to three. The six non-zero coeffi- 
cients of this polynomial can be computed from T only. Therefore evaluation of this 
term only requires the evaluation of this polynomial. 
For this functional the error bounds 

\V{T) - V{T)\ < cih'^U\V{T)\ 

hold, with analogous bounds for the volume of the whole triangulation. Here h is 
the maximum of the diameters of T and T, U — inax{\u(pi)r(pi)\ : i = 0,1,2} and 
ci = ci{g, dg, d^g) depends on the metric with ci < Cr~^ in view of definition 12. 61 

The discrete volume functional is again as differentiable with respect to the nodal 
values of u as the metric g. Explicit formulas for these derivatives can be computed 
from in the same way as for the area functional. 

The discrete volume and area functionals are therefore constructed to give the 
right answer in the Euclidean case since the asymptotically flat situation suggests 
that this gives a good guess. Indeed, the error bounds given above improve rapidly 
for growing radius. Their values and derivatives can be computed per triangle using 
a fixed procedure, that is fixed time. Thus evaluation of these two functionals takes 
time proportional to the number of triangles. 

Evaluation of the functionals and their derivatives requires the evaluation of the 
metric and its derivatives at the requested points. Therefore, if the metric is given on 
a grid it is necessary to interpolate these values. 

Note that the gradient of the area functional is a discrete, weak analogon to 
the mean curvature. This can be seen from lemma and the fact that deforming 
graphu over a triangulation by increasing Up at one point p corresponds to a variation 
of graphw with variation vector field (ppVp. 
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3.3. Optimization techniques 

Now we have a problem of the foUowing form: 

minimize A{u) u G 5 ~ 

s.t. V{u) - Vo 

with A{u) the area of graphu over a fixed reference triangulation and V{u) the discrete 
oriented volume between graphw and the reference triangulation as described above. 

Naturally one would consider a Lagrange method, that is computing the critical 
points of the Lagrangian function 

L : 5 X R ^ R : (m, A) i-> A{u) - A (v{u) - Vo) ■ 

However, since a local minimum of Q does not correspond to a minimum of L but 
merely to a critical point of L, minimization methods are not applicable. Thus solving 
the critical point equation directly with a Newton method is the only alternative in 
view, but this is not desirable, since it would involve second derivatives of the metric 
9- 

For numerical optimization however, there is a better method, namely the 
augmented Lagrangian method. This considers the following penalized Lagrangian 
function 

Lp{u, A) - A{u) - A (V{u) -Vo)+^ {V{u) - Vo) ' 

with a penalty parameter p. The following theorem can be proved using standard 
calculus methods on the submanifold generated by the constraint. However, a more 
elementary proof can be found in |15[ Theorem 12.2.1]. 

Theorem 3.1 If A and V are , u* is a solution to A* is the exact value of the 
Lagrange parameter, and the Hessian of A at u* is positive in directions perpendicular 
to gradV , then u* is a critical point of Lp{-,\*) and there exists po .such that for all 
p > po the Hessian of Lp{-, X*) is positive definite, that is, u* is a strict local minimum 

0fLp{;X*). 

Note that in the analytic case, theorem 12.71 implies that the conditions of this 
theorem hold. A CMC-surface is therefore a local minimum for the analytic penalized 
Lagrangian for suitable penalty and Lagrange parameters. For the discrete case we 
can still assert the regularity assumption, but we unfortunately do not know about 
stability. 

To solve (@J) we minimize the augmented Lagrangian. To find the Lagrange 
parameter to the desired volume the following algorithm is used: 
p <— some value > po 
Ao ^ good initial guess 
fc ^ 
repeat 

Minimize ip(-, Afe) to get approximative solution Uk 

k^k + i 

until |A - A*| < e 

If the conditions of theorem 13 . II hold, then this algorithm converges locally in (A^) to 
A* and then the {uk) converge to u* . Convergence improves for p — > oo. For a proof 
of this fact cf . JS| ■ 
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A "good initial guess" for A can be obtained from the Euclidean situation. If a 
surface is considered that has a radius of approximately r then 2/r, the Euclidean 
mean curvature of a sphere of radius r, is a good choice, but in practice, starting with 
A = also works well. 

Each pair [uk, Xk) obtained by one step of the above algorithm corresponds to a 
discrete CMC-surface and its discrete mean curvature, of course not with V{uk) = Vq, 
but rather V{uk) = V{uk). This explains why not many steps have to be performed 
when one just wants to find CMC-surfaces together with their mean curvature and 
area, but not necessarily a particular enclosed volume. 

The parameter p should neither be chosen too small nor too big, since a value 
too small gives that the minimization method diverges and a value too big gives a 
bad condition of the minimization problem. An estimate for pQ can be obtained by 
estimates for the Jacobi operator from Huisken-Yau . These estimates indicate that 
Po decays like r~^, where r is the radius of the surface considered. But the condition 
of the problem is not affected if p is chosen to decay like since that is the scaling 
of the other terms of the Hessian matrix of Lp, which controls the condition of the 
problem. Therefore the exact choice of p is not very critical, but decreasing p should 
be attempted. 

The penalized Lagrangian method is used to transform the constrained problem 
into an unconstrained minimization. The method we chose to numerically solve this 
is the conjugate gradient method. Although the conditions of theorem 13.11 imply 
local convergence of such a method, the rate of convergence depends on the basis 
chosen for iS. The nodal basis is not a very good choice, since the number of steps to 
make increases proportionally to the number of points of the triangulation and would 
therefore give an algorithm with quadratic time complexity. 

The similarity of the problem to the solution of linear elliptic PDEs suggests that 
one try methods from this field that have proven to give good convergence rates. The 
choice made here is to use hierarchical bases that give the CG method multigrid like 
speed while being very simple to implement. Using this preconditioner, an overall time 
complexity of ©(TV log TV) for linear problems is achieved with N being the number of 
vertices of the triangulation. Section H31 describes some examples and shows that the 
speedup is significant, and reduces complexity even in the nonlinear case. 

Optimization methods are presented in Fletcher |17l I15| . For details on 
hierarchical bases see Bank, Dupont and Yserentant ^H] and Yserentant |19l 120! • 

3.4- Computation of single surfaces 

When a single surface has to be computed we can take the full advantages of the 
hierarchical finite element representation, and use a cascading technique for iteration 
as proposed by Bornemann and Deufihard [23 • We start with a given coarse 
triangulation and iterate to get a first coarse approximation to the surface. Then 
we refine this triangulation by dividing each triangle into four new triangles. The 
surface obtained from the coarse grid iteration then can be used as initial value to 
the fine grid iteration. This procedure can be continued until the desired resolution is 
reached. This cascading technique gives an immense speedup. 
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Figure 2. A translated triangulation has reduced resolution. This is the final 
iterate of an area minimization in Schwarzschild, started from a translated sphere. 

3.5. Computation of foliations 

When the method is used to compute whole famihes of constant mean curvature 
surfaces close to Euclidean spheres it is possible to use the previously computed surface 
as initial data for the next iteration. To do that, one has to produce a sequence of 
appropriate reference surfaces and volume conditions that give a reasonable sequence 
of surfaces. 

The following is suggested in a situation when one knows one discrete constant 
mean curvature surface T that is centered. In this case, the radial direction on this 
surface can be used as the variation direction of the vertices. 

To produce a new surface at distance r further out (in), T should be replaced by 
the graph of the constant function r (— r) over T. This graph then will be the new 
reference surface T. 

The constrained minimization problem to consider now is to minimize A while 
keeping V" = 0. This leads to an algorithm where not much volume is inbetween the 
reference surface and the unknown constant mean curvature surface, which is desirable 
since this increases the accuracy of the discrete volume functional. 

The problem of finding the first constant mean curvature surface can sometimes 
be solved by simply minimizing area without the constraint. This will give a discrete 
minimal surface, which can serve as starting surface. However, not all manifolds have 
a single spherical minimal surface that can be used. 

Another method to get a starting surface is to start with a Euclidean sphere of 
big radius, use its normal direction as variation direction and solve the constrained 
problem. However, if a significant translation occurs, the solver should be restarted 
using a coordinate sphere around the computed center as initial surface, since the 
distortion in the triangulation reduces the resolution, for illustration see figure El 
Alternatively one could use adaptive mesh refinement here. 

3.6. Ceneralizing the Algorithm 

The algorithm described above is based on the variational structure of the equation 
H — const . 

In view of numerical relativity it would be very interesting to find surfaces S that 
solve the equations 

H±P = const (4) 
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with P — tiY, K — tvM K — Klv^v), v the unit normal to S, and {M,g,K) are an 
initial data set for the Einstein equations. 

However these equations are not the Euler-Lagrange equations for any functional 
considering surfaces in M . This is due to the fact that the term P depends on the 
normal of the surface in consideration. In contrary the equation 

H + Hf) = const 

where Hq : M ^ R does not depend on the normal of the surface is the Euler-Lagrange 
equation for a constrained minimization problem, in the sense that it is satisfied on 
the boundary dft of a set il C M minimizing 

j{n) = \dn\ + [ Hodfi n<zM 

Jn 

subject to / dfj, = Vo = const . 
Jn 

For a given initial surface Sg with normal vector field vq we may therefore fix a 
function Pq — tr K — K^vq, vq) and extend it to a neighborhood of the initial surface 
by the condition that it does not change along the radial directions considered in 13. II 
Then we solve the equation H ± Pq — const using the procedure described before. 
An iteration of this strategy will produce a sequence of surfaces and normals that 
might converge. If they converge, then the limit is a solution of the discrete version 
of equation^ 

The author currently works on implementing, testing and examining this 
approach, but came to the conclusion that this preliminary outline for extending this 
algorithm might be interesting for applications in numerical relativity. 

4. Numerical Examples 

In this section we present three examples based on metrics of the Brill-Lindquist-Type. 
This is a conformally flat metric on of the form 

This metric represents a spacelike slice containing N Schwarzschild-Type singularities 
at the points Xk of mass . The metric and its derivatives were evaluated analytically 
at every requested point. 

The triangulations used for the computations in this section are based on the 
octahedron, that is the triangulation with vertices ie^, i = 1,2,3, and regular 
refinement of it. To regularly refine a triangulation, every triangle is divided into 
four subtriangles by introducing the midpoint of its edges as new vertices to the 
triangulation. If these new points are projected to the sphere one obtains the surfaces 
shown in figure 13 Rescaled and translated versions of these triangulations in different 
refinement states will serve as starting surfaces. 

All pictures of surfaces shown in this section were created using geomview |22| . 

4.I. Schwarzs child Solutions 

Here iV = 1, mi = 1 and xi — 0. This case was included to test the convergence of 
the method. To compare the numerical results with the analytically known results. 
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Figure 3. Several refinement states of tlie octahedron. 
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Figure 4. Exact and com- 
puted mean curvature com- 
pared for different radius. 
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Figure 5. The difference of 
the Hawking mass for differ- 
ent translated initial spheres 
and the expected value 1. 



introduce the following radius function for a surface S 



= |S| 



\x\d^i 



where | • | denotes Euclidean distance to the origin. Using the center of gravity rule 
for numeric integration, this radius can be computed approximately, which gives for 
each triangulation T a discrete radius r(T). 

Starting the method with the minimal surface and computing outward, one 
obtains a family of triangulations , and a family of Lagrange parameters which 
can be interpreted as the constant mean curvature of these triangulations. The exact 
mean curvature of a sphere of radius r in the spatial Schwarzschild metric is plotted 
as continuous line in figure 01 whereas the Lagrange parameters versus the numerical 
radius of each triangulation is plotted as mark. This data results from a resolution 6 
triangulation. 

Since a good approximation to the mean curvature is known, the Hawking mass 
on constant mean curvature surfaces 

|y|i/2 

"^h(S)- 77^77 (16^-15^1^') 



(16^)3/2 ^ 

can easily be computed. Figure |S1 shows the absolute value of the difference of the 
Hawking mass of each surface and the expected value 1. The initial surfaces were not 
centered here but translated by d = 0, 0.1, 0.2, 0.3, 0.4, resolution is again 6. 

To test the ability of finding the center, the method was started with Euclidean 
spheres of radius 20 and center {d, 0, 0) for different values of d. Then a CMC surface 
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Table 1. Distance of the computed center to the origin for translated starting 
surfaces of radius 20 off center by distance d and the associated convergence factor 
/ for d = 5. 

Resolution d=l d = 2 d = 3 d = A d = 5 f d = 5 

4 0.01791 0.03521 0.05394 0.07245 0.09161 4/5 3.93 

5 0.00455 0.00903 0.01382 0.01869 0.02329 5/6 3.95 

6 0.00112 0.00230 0.00343 0.00472 0.00589 6/7 3.98 

7 0.00030 0.00060 0.00086 0.00115 0.00148 





40 60 
Iteration 



Figure 6. A rough initial surface on the left and the the maximum of the absolute 
deviation of r from 1/2 during an iteration for this surface. For comparison the 
data for an iteration with smooth data is shown. 



of equal volume was computed using the value 0.0005 as penalty parameter. The 
Euclidean distance of the center of this surface to the origin is shown in table ^ 
for different resolutions and different values of d. The resolution number is the 
number of refinement steps applied to each triangle. Resolution 3,4,5,6,7 correspond 
to 258,1026,4098,16386,65538 vertices. This table shows nearly the expected second 
order convergence, which corresponds to a factor of 4. 

A direct elliptic method fails to converge when started with a very rough surface 
near the horizon, as is reported by Thornburg To test this aspect of our method, 
we construct rough data by starting with a sphere of radius 1/4 and refining up to 
level 5. The vertices that are created by the last refinement step are moved outward 
to a radius of 3/4 such that the starting surface oscillates around the minimal surface. 
In comparison to spheres of radius 1/4 and 3/4 as starting surfaces, the iteration takes 
significantly longer, about three times, but still converges. The diagnostic quantity 
E := max^^Y vertices I'' ~ '^■''1 plotted in figure|n|for spheres of radius 1/4 and 3/4 
and the rough surface as initial surface. The iteration stopped, when the relative 
change in E was less than 1/1000. 

4-2. Two singularities 

Here we examine a metric with N = 2, — 1 and Xk = ±c?ei. This space contains 
two minimal surfaces, each one enclosing one of the singularities. 

The first test was to find the critical value d* , such that the two singularities are 
enclosed by common a minimal surface for d < d* but not for d > d*. Refinement 
up to level 7 indicates that \d* ~ 0.766362| < 0.000002, which deviates by 0.02% from 
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Figure 7. The error the 
Hawking mass for two singu- 
larities of mass 1. 



Figure 8. The mean 
quadratic deviation of the ra- 
dius function for two singu- 
larities of mass 1. 



Table 2. Area, mean curvature H, Hawking mass m and associated convergence 
factors for the surfaces computed form a coordinate sphere with radius 5 in the 
d = 1-metric for two singularities. 

Resolution Area H m /Area /h /™ 

3 646.68121 0.1860585 1.9893651 3.97 3.88 3.80 



4- 5 

5 650.11737 0.1856342 1.9934694 $^ 4.00 3.97 3.95 

5— 6 



4 649.42623 0.1857211 1.9926138 

5 650.11737 0.1856342 1.9934694 

6 650.29032 0.1856123 1.9936860 



the value 0.76619745 reported by Thornburg in j^. The area of this minimal sm^facc 
is computed as 196.41579, whereas Thornburg reported a value of 196.407951 which 
deviates by 0.004%. We noticed however, that with increasing resolution the value 
of d* and the area of the minimal surface decreased with our procedure, therefore 
for larger resolution the values reported here might become closer to the values of 
Thornburg. 

For d > d* , computation cannot start at a minimal surface, since this does not 
exist. That is why for this problem the algorithm is started from a Euclidean sphere 
with big radius centered at 0. In each step the radius has been reduced approximately 
by 1/10. In figure El every tenth surface for d = 5 is drawn. The two small surfaces 
are the minimal surfaces enclosing the singularities. 

It is clear that for growing radius, the Hawking mass of the surfaces approaches 
the ADM-mass, which in this case is the total mass of the two singularities, namely 
2. The difference of of the Hawking mass and the expected value is plotted in figured 
for d = 1,2,3,4,5. 

The computed area, mean curvature and Hawking mass for d = 1 of a surface 
with approximate radius 5 for different resolutions and the associated convergence 
factor is shown in tabled 

Figure d shows the mean quadratic deviation of the Euclidean radius function 
from the expected radius. This figure shows that the surfaces become rounder with 
growing radius as predicted by the theory. 
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Figure 9. The foliation for two singularities of mass 1 and distance 5 to the 
origin. The innermost surface is the last surface computed by the method, for 
smaller surfaces the algorithm seems to fail converging. 

4-3. Three singularities 

Here we study a very symmetric Brill-Lindquist metric with three singularities at 
Xi = 0, X2/3 = ±ei and masses mi = 2, TO2 = — 2\/2- This comes from the 
stereographic projection of a metric conformal to the standard metric on C R"*. 
The conformal factor on the sphere is the sum of the Greens functions to the conformal 
Laplacian at the points ±(1, 0, 0, 0) and ±(0, 0, 0, 1), rescaled such that in the standard 
stereographic projection of S'^ \ (0, 0, 0, 1), the resulting metric has the Brill-Lindquist 
form. 

An analogous picture is shown in figure E| where it is clarified that for each 
asymptotically flat end an outermost minimal surface exists. Additionally there are 
two minimal surfaces resulting from the mirror symmetry along the planes, that divide 
the picture in equal halves. As it turns out, these additional minimal surfaces are stable 
and can be found by minimization. 

The intersection of these minimal surfaces with the x^x^-plane is drawn in 
figure El in the stereographic projection. The smaller surfaces enclosing one 
singularity and the big surface enclosing all three singularities correspond to the 
outermost minimal surfaces of the asymptotically flat ends. The surfaces enclosing 
two singularities correspond to the surfaces arising from the symmetry. The area 
of the outermost minimal surfaces is given in figure |3| The surface named (R)ight 
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Figure 10. 5^ with four points blown to infinity by conformal rescaling, 
illustrated by on the left and in the stcreographic projection with minimal 
surfaces on the right. 



Table 3. The area of two different minimal surfaces and the convergence of its 
difference. 



Resolution 


R 


M 


O 


f 


R-M 


R-O 


3 


2966.95 


2970.59 


2970.84 


3/4 


4.10 


4.09 


4 


2951.64 


2952.53 


2952.59 


4/5 


4.02 


4.02 


5 


2947.79 


2948.01 


2948.03 


5/6 


4.00 


4.00 


6 


2946.83 


2946.88 


2946.89 


6/7 


4.00 


4.00 


7 


2946.59 


2946.60 


2946.60 









encloses the singularity at (1, 0, 0), the surface called (M)iddle encloses the singularity 
at the origin, and the surface called (O)uter encloses all three singularities. These 
areas should be the same due to symmetry. Therefore we compute the convergence 
factors of the differences to zero and obtain, as also shown in figure perfect second 
order convergence. 

This symmetry enables us to compute these surfaces exactly. In R** the surfaces 
are the intersections of with the planes = and project to the surfaces given 
by the quadratic equation 

Q{x) := \xf ± 2a;^ - 1 = 

In figure ITTI we plot maXg^jj vertices I '3(^)1 ^'^^ iteration started on a sphere 
with radius 1 centered at (0.5,0,0) for different resolutions. We again find quadratic 
convergence for increasing resolution. 

4- 4- Kerr Solutions 

For this test, a {t = const}-slice of the Kerr metric in Kerr-Schild coordinates is used. 
The associated 3-metric is given by 

with 

MR / Rx^ + ax^ Rx^ - ax^ x^ \ 

^ " i?2 + a2 cos2 Q i?2 + a2 ' _^ ' j 
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Figure 11. The maximum value of \Q(x)\ for iterations in different resolutions. 
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Figure 13. The mean curvature of the CMC-surfaces in Kerr for a = 0.3, 0.5, 0.7 
and 0.9. 



where 

i?2 = i(r2 - a'^) + ^i(r2 - a^y + 02(^.3)2 cos^ 6* = 

and jxp is the EucHdean distance. This metric has an apparent horizon at 

{R — R+}, R+ = M + V — if a < 1. Unfortunately in this shcing the horizon is 
not a CMC-surface in contrast to the Boyer-Lindquist {t — const} shces, and cannot 
be found with the method presented here. 

We therefore start again at a large radius r — 20 with a centered sphere and 
compute inwards reducing the radius by 0.1 in each step. Some of the inner surfaces 
for different values of a are shown in figure^l The last surface shown is the last surface 
computed by the method. The mean curvature of some more surfaces is plotted for 
these values of a in figure [T51 while the mean quadratic deviation of the radius function 
from constant radius and the hawking mass is shown in figure [T^ for the same values 
of a. This figure displays the remarkable fact that the hawking mass of the CMC- 
surfaces is nearly independent of a, while not constant. The same shape of the graph 
even holds for a = 0, where the CMC-surfaces are perfectly round coordinate spheres. 
Figures El and E| show the data of 190 surfaces that are nearly equidistant with 
distance 0.1, while in figure [T^ onlv 15 of the innermost surfaces with distance 0.2 
are shown. The apparent horizon is shown in none of these pictures, since it is not a 
CMC-surface. In the above plots we use the geometric area radius Vg — ■y/|E|/47r on 
the horizontal axes. 

4-5. Performance 

This section reviews the performance of the algorithm. All times reported here have 
been measured on a personal computer with a single Athlon XP 2000-)- CPU. The 
first performance test is intended for comparing this method to others. The metric 
used is the unit mass Schwarzschild metric gs = (1 + ^)^<5- The aim was to find the 
horizon, so area minimization without volume constraint was attempted. The initial 
surface is taken as discrete sphere of radius 1 with center (0.3,0,0). We use a low 
resolution of 3 to start and initially optimize. Then we refine to the next resolution, 
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Table 4. Elapsed Time and Iterations for the Schwarzschild test. 



Resolution 


Number 


Iterations 


Total 


final value 




of points 


on finest Level 


Time [s] 


of E 


3 


258 


25 


0.44 


0.014472349 


4 


1025 


22 


3.22 


0.003732145 


5 


4098 


26 


16 


0.000945821 


6 


16386 


21 


58 


0.000236132 


7 


65538 


39 


374 


0.000060190 


8 


262146 


42 


1726 


0.000015141 



take the surface obtained from the coarse grid iteration as initial values for the fine 
grid optimization, and iterate this process until the desired resolution is obtained. 

Since the horizon is a Euclidean sphere of radius 1/2 we use E :— 
max^jj vertices l'^"-'^/^! diagnostic parameter. Iteration is stopped when E's relative 
change is less than 1/1000 for more than four iterations. Table0]shows the iterations 
elapsed for each refinement level and the total time elapsed up to reaching the result, 
including the coarse grid iterations. Note that the final values for E display nearly 
perfect second order convergence. Figure El shows the value of E for each iteration of 
the algorithm with maximal refinement level 8 starting with refinement level 3. The 
levels of the cascade correspond to the convergent regime of the algorithm for the 
current refinement level, while the edges occur when refinement has been done. 

To test the speedup of the hierarchical basis transformation, we go back to Brill- 
Lindquist data with two singularities and d = 1. We start with a large coordinate 
sphere of radius 10, translated by 1 and test the time used to center this sphere to 
a constant mean curvature surface using the nodal and hierarchical basis. We again 
used the cascading technique with initial refinement level 3. As diagnostic parameter 
we use the the ^^-gradient norm times two to the power of the number of refinement 
levels, which is nearly independent of the refinement level. The stopping condition 
was to reduce this value by a factor of 2000. We chose to do only one minimization 
and therefore started with the previously computed Lagrange parameter A — 0.13531. 
The total time elapsed and the number of iterations in the highest resolution is plotted 
in figure 1161 We see a substantial improvement in the number of iterations involved 
and the elapsed time. 
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Figure 15. The maximum of the absolute value of the radius function minus 0.5 
monitored during the iteration of a three to eight refinement steps run. 
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Figure 16. Total time (left) and the number of iterations for the highest 
refinement level (right) plotted versus the number of refinement levels. 



4-6. Comparison to other methods 

Comparison of the speed of this method to other methods available is very difficult 
at this stage. This is mainly due to the following three issues. At first there are no 
standardized tests nor standardized platforms on which to test nor agreements with 
what accuracy to stop, which are general problems. Second, the figures presented by 
Schnetter |S| and Thornburg do not refer to the computation of CMC-Surfaces 
but to the computation of apparent horizons, which is a different equation in general 
and especially in the cases considered by them. And third, as stated before, our 
method uses an explicit evaluation of the metric, whereas Schnetter and Thornburg 
use interpolation of the metric from a given grid discretizing the metric on the ambient 
three-manifold. As evaluation of the metric takes more than half of the execution time 
in our case, it is hard to judge whether we compare the underlying algorithms or just 
the different methods of evaluating the metric. 

However, in table |S1 we display some figures taken from Thornburg to show in 
comparison to tabled] at least that the execution time of our algorithm is in the same 
order of magnitude as the times given by Thornburg. Thornburg tries to find the 
horizon in a boosted slice of Kerr data with angular momentum a = 0.8 and velocity 
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Table 5. Values reported by Thornburg for horizon finding in a boosted Kerr 
time slice, with different resolutions. 



Number of points 


Total Time[s] 


533 


2.0 


1121 


4.2 


2945 


13 


7905 


43 


25025 


220 



V = 0.8. The numbers refer to execution time on a dual Intel Pentium IV 1.7 GHz 
which executes his algorithm on one processor, where in contrast we worked on a single 
AMD Athlon XP 2000+ . 



5. Conclusion 



We have presented an algorithm to compute constant mean curvature surfaces based 
on a finite clement discretization. The implementation displays convergence, accuracy 
and speed rivaling the methods based on finite differencing used in production runs. 
But still the implementation by the author is academic in the sense that not much 
effort has been spent on optimization and fitting the code into standard environments 
such as Cactus for example. 

However, creating a fully functional apparent horizon finder based on finite 
element discretization, i.e. the extention of this algorithm to not necessarily time 
symmetric initial data, is still a long way to go, but the author's opinion is that the 
results presented here are very encouraging. 
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